Thermodynamics of tf-wave superconductors in a magnetic field 



O 

o 
o 

(N 
> 

o 



o 
o 

c3 



-a 
c 

o 
o 



On 
O 



O 
O 
-i— > 



i 

C 

O 

o 



I. Vekhter 1 !, P. J. Hirschfcld 2 , E. J. Nicol 3 
1 Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545 
2 Department of Physics, University of Florida, Gainesville, FL 32611 
3 Department of Physics, University of Guelph, Guelph, Ontario NIG 2W1, Canada 

(February 1, 2008) 

We investigate the thermodynamic properties of a two-dimensional d-wave superconductor in the 
vortex state using a semiclassical approach, and argue that such an approach is valid for the analysis 
of the experimental data on high-temperature superconductors. We develop a formalism where the 
spatial average of a physical quantity is written as an integral over the probability density of the 
Doppler shift, and evaluate this probability density for several model cases. The approach is then 
used to analyse the behavior of the specific heat and the NMR spin-lattice relaxation rate in a 
magnetic field. We compare our results with the experimental measurements, and explain the origin 
of the discrepancy between the results from different groups. We also address the observability of 
the recently predicted fourfold oscillations of the specific heat for the magnetic field parallel to the 
copper oxide planes. We consider both the orbital and the Zeeman effects, and conclude that at 
experimentally relevant temperatures Zeeman splitting does not appreciably reduce the anisotropy, 
although it does change the field dependence of the anisotropic specific heat. We predict a scaling 
law for the non-exponentially decaying NMR magnetization, and discuss different approaches to the 
effective relaxation rate. 



I. INTRODUCTION 

Despite significant recent advances we still lack a com- 
plete understanding of the physics of low energy excita- 
tions in the vortex state of unconventional superconduc- 
tors. High-temperature superconductors (HTSCs) are an 
example of a system where theoretical predictions can be 
checked against a large bodof experimental evidence. In 
zero field these materials have a d-wave superconduct- 
ing energy gap, with nodes along the diagonals of the 
Brillouin zone, and consequently a finite density of low 
energy excitations.El Moreover, it is believed that at tem- 
peratures low compared to the transition temperature, 
T -C T c , these excitations are reasonably well described 
by the Landau quasiparticles, even though such an ap- 
proach fails in these materials at higher energies. A vari- 
ety of experimentally measured quantities such as the 
electronic j-specific heatiSu, effective penetration depth 
from /iSRO, [Soia-lattice relaxation rateLrB, and thermal 
conductivityEnLil are available to test the predictions of 
theories. 

In this work we discuss the influence of the magnetic 
field on the thermodynamic quantities in the vortex state 
of the unconventional superconductor, and, in particular, 
address the question of how these properties depend on 
the structure of the vortex state. We concentrate on the 
behavior of the density of states and the electronic spe- 
cific heat and the spin-lattice relaxation of the NMR mag- 
netization. There exist several theoretical approaches to 
the analysis of the thermodynamic quantities in the vor- 
tex state of unconventional supepennductors. We employ 
here a semiclassical approaclJljO, which has been suc- 
cessful in describing tha-field dependence of a variety of 
the physical quantitiesll3ilJ. It is an approximate de- 
scription, and in the next section we discuss the region 



of its validity and the grounds for our belief that it is 
applicable to the present problem. Section III introduces 



the basic model of the nodal quasiparticles and the idea 
that the physical quantities in the vortex state can be 
obtained by calculating the spatial average of their lo- 
cal values, computed with the help of the semiclassical 
approach. 

Until now such spatial averages have only been done 
analytically in an oversimplified model of a single 
vortexEjO In this paper we introduce a generalization 
of these approaches by rewriting the spatial average as 
an average over the probability density of the Doppler 
shift of the quasiparticle energy in the presence of the su- 
perflow. Restating the problem in this language enables 
us to introduce several model distributions of the prob- 
ability density, discussed in Section [V, and investigate 
how the physical quantities obtained within the semi- 
classical framework depend on these distributions and on 
the structure of the vortex state. We obtain the energy 
and field dependence of the density of states for the ge- 
ometries with the magnetic field applied both normal to 
the superconducting planes and in the plane in Section 
0. This density of states is used to analyze the behav- 
ior of the electronic specific heat in HTSCs. We obtain 
the energy scales relevant to the high-temperature super- 
conductors in the vortex state, and suggest a resolution 
to the origin of the disagreement between different ex- 
perimental groups regarding the magnitude of the field- 
dependent term in the specific heat and the form of the 
scaling function; this is the content of Section VI. In the 



same Section we address the question of the observabil- 
ity of the oscillations in the specific heat for the magnetic 
field applied in the superconducting plane as a function 
of the angle between the field and the nodal ditecjions. 
These oscillations have been recently predictecESlIa, but 
so far have not been observed. Part of the difficulty may 
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stem from the smallness of the in-plane Doppler energy 
scale, as inferred from the experimental measurements; 
it has recently been argued that the Zeemaa-splitting 
reduces the observed oscillations significantlyjij and we 
investigat e its effect in detail in this work. 

Section VII is devoted to the effect of the nonuniform 



density of states on the spin-lattice relaxation time. This 
non-uniformity leads to a non-exponential decay of the 
magnetization .and to a field dependence of the effective 
relaxation rateEj; here we show that the effective relax- 
ation rate depends on the structure of the vortex state, 
and obtain an approximate form for it. We also predict 
a scaling law for the magnetization decay which can be 
checked directly. The effect of impuriti es on the density 
of states is briefly addressed in Section VIII . We expect 
these effects to be very important for the discussion of 
the transport properties which we defer to a later pub- 
lication. Finally we summarize our findings and discuss 
some open questions. 



II. SEMICLASSICAL APPROXIMATION 

HTSCs are extreme type II superconductors (the ra- 
tio of the London penetration depth Al to the coherence 
length £o is large, Xl/£,o ~ 100), and are in the mixed 
state over the range of applied fields, H, from a few hun- 
dred Gauss to well in excess of 50 Tesla in YBa2Cu307_,5 
(YBCO) and Bi 2 Sr2GaCu 2 8+5 (Bi-2212) near optimal 
doping. In the mixed states the magnetic field penetrates 
the bulk of the superconductor in the form of vortices, 
which consist of the cores, where the superconducting 
order parameter is suppressed, and circulating supercur- 
rents around them. The vortex core-size is of the order 
of the coherence length, £o ~ 15AEj ! til, while the average 
intervortex distance can be estimated by imposing the 
requirement of one flux quantum $0 = hc/2e per vor- 
tex, or d/2 = R = yj $0/71- £?, where B is the internal 
field. At typical experimentally accessible fields (1-20 T) 
\l 2> d ^> £oi the magnetization due to the vortex lat- 
tice is small, and the internal field can be replaced by 
the applied field, H, so that dy/H ~ 500A-T 1 / 2 . The 
actual distance differs from the average value d by a nu- 
merical factor of the order of unity, which depends on 
the structure of the vortex state; the vortices in HTSCs 
may form a regular lattice (as they do in YBCO and in 
Bi-2212 at low fieldsua) orJae moderately disordered 
(as in Bi-2212 at higher fieldsO). 

In the experimentally relevant fields, B c \ <C H <C H C 2, 
where H cl (H C 2) is the lower (upper) critical field, there 
may exist two types of low energy excitations. First, as 
in conventional, s-wave materials with an isotropic gap, 
there may be a branch of law.-energy fermionic excitations 
bound to the vortex coresEj Experimental evidence sug- 
gests, however, that there is at most .one such state in the 
vortex cores of YBCO and Bi-2212t3cil, and therefore the 
properties of the mixed state of d-wave superconductors 



are dominated by the "extended" quasiparticle states in 
the bulk. These states are formed when quasiparticles 
with momenta close to the position of the nodes of the 
gap, Ak, in the momentum space (and therefore with a 
small gap) interact with the supercurrents in the vortex 
state. Most of the theoretical work has therefore explored 
the properties of these states. 

Very significant progress has been made by utilizing 
the semiclassical approach, which treats the momentum 
and the position of the quasiparticle as commuting vari- 
ables. It is valid when the wave function of a quasipar- 
ticle can be replaced by its envelope on the length scales 
exceeding the coherence length, i.e. when fc/£ 3> 1, 
where kt is the inverse Fermi wavelength. In that method 
the effect of the supercurrents is accounted for by in- 
troducing-a-pPoopkr shift into the quasiparticle energy 
spectrum,Ea€3'BB E'(k, r) = E(k) + e(k, r). Here E(k) 
is the energy of a quasiparticle with momentum k in the 
absence of the field measured with respect to the chemi- 
cal potential. In a two-dimensional <i-wave superconduc- 
tor this spectrum is conical (massless anisotropic Dirac 

spectrum): E(k) w ± x Jv^k? L + i> A fcjj, where the Fermi 

velocity Vf is associated with the dispersion of the quasi- 
particles in the direction normal to the Fermi surface 
(component k± of the momentum), while ua ~ Ao/fe/ 
is the slope of the gap at the node associated with the 
dispersion of the quasiparticles along the Fermi surface 
(fell). The Doppler shift, e(k,r) —v s ■ k depends on the 
quasiparticle momentum and the local value of the su- 
pervelocity, v s (r). This shift in the energy is an exact 
result for a uniform supercurrent,EZI where it reflects the 
pairing of the electrons with a finite center of mass mo- 
mentum. In the simplest picture such an approach re- 
mains valid for a non-uniform current for as long as the 
spatial variations of v s are slow on the scale of the spa- 
tial extent of the Cooper pair, £o- In superconductors 
with nodes in the energy gap, the Doppler shift may ex- 
ceed the local (in the momentum space) gap, and leads to 
an increase in the density of the unpaired quasiparticles: 
even at T — for some positive energies E the shifted 
energy E' is negative so that the corresponding states 
become occupied. In the context of d-wav.e-superconduc- 
tors this was emphasized by Yip and SaulsEa, who investi- 
gated the effect of the screening-currents in the Meissner 
state on the superfluid density.Ea These currents vary on 
the scale of the penetration depth Al ^> £0, so that the 
Doppler shift description is appropriate, and result in a 
linear dependence of the effective penetration depth .on 
the applied field for certain experimental geometriesEa 
However so far the predicted dependence has not been 
confirmed experimentallyES 

Similar physics is at play in the dilute vortex limit. At 
distances small compared to the penetration depth, su- 
percurrents around an isolated vortex are inversely pro- 
portional to the distance from the center of the vortex, r; 
for £0 <C r <C Al the supervelocity field is |v s | = h/2mr, 
where m is the quasiparticle mass. Consequently, the re- 
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quircment of the slowness of the variation of v s , which 
can be written as |Vv,,|£o << v s (two particles compris- 
ing the Cooper pair "see" the same velocity), is satisfied 
at r 3> Co i justifying the use of the semiclassical approach 
outside of the core and therefore for the analysis of the 
extended quasiparticle excitations at energies low com- 
pared to the gap maximum. Such an analysis for a single 
vortex in <i-wave superconductors was first carried out 
in a landmark paper by Volovikllj, who showed that the 
density of the extended quasiparticles at low tempera- 
ture, T, varies s,s y/W. This result was confirmed 
by Moler et alB and subsequently by other groups^ 
from the measurements of the electronic specific heat in 
an applied magnetic field. Numerical studies of the tight - 
binding-|inodel are also in qualitative agreement with this 
result E3 Moreover Volovik has shown that the density of 
the extended quasiparticles dominates that of the states 
bound to the vortex-core even if the latter set is treated as 
a quasi-continuun£3 (as it would be in a superconductor 
with a long coherence length) , providing further theoret- 
ical support for neglecting the core states in the analysis 
of the properties of the vortex state in unconventional 
superconductors. 

The semiclassical approach was incorporated into the 
Green's function formalism by Kiibert and HirschfeldEJ, 
and used in that form to analyze thermodynamic and 
transport^rxupcrties of the high-T c cuprates in the vor- 
tex state. Edil3 In particular, accounting for the impurity 
scattering in this framework has significantly improved 
the agreement between the theoxy and the measurements 
of the electronic specific heatEj, and the field depen- 
dence of the low-temperature thermal conductivity^ is 
in qualitative agreement with the results of a semiclas- 
sical calculations. In the semiclassical approach the ef- 
fect of the magnetic field is contained in the new en- 
ergy scale associated with the Doppler shift, Eh — vj/d, 
and the behavior of the physical properties is determined 
by the competition between this energy, the tempera- 
ture, and the impurity scattering rate. Phot popi ission 
measurements on high-T c compounds suggestt3Eil that 
v f ~ 1.5 - 2.5 x 10 7 cm/s leading to E H ~ 
K-T- 1 / 2 . 

For a long time, understanding of the low-energy exci- 
tations in the vortex state beyond this semiclassical pic- 
ture has proved elusive. The difficulties stem in part from 
the need to treat on equal footing the applied magnetic 
field and the superconducting currents (semiclassical ap- 
proach treats the supercurrents classically). Attempts 
have been made to take as a starting point the Lan- 
dau quantization of the quasiparticle states^ and include 
the effects of supercurrents perturbativelyJ£3 : EJ however, 
since the supervelocity field is long-ranged and singular at 
the position of each vortex the Landau levels-are strongly 
mixed, making a detailed analysis difficult.Efl 

The most significant progress has ibeen made in a 
recent work by Franz and TesanovicJ^j who have in- 
troduced a gauge transformation which takes into ac- 



count both the supercurrent distribution and the mag- 
netic field. In their approach the problem is mapped 
onto that of nodal Dirac fermions in an effective zero 
average magnetic field interacting with effective scalar 
and vector potentials which are periodic in the uuijt cell 
of the vortex lattice. Both Franz and TesanovicLj and 
Marinelli et alx3 have studied the band structure of the 
nodal quasiparticles for perfectly periodic vortex lattices 
for various values of the anisotropy of the Dirac spec- 
trum, an — Vf/vA- 

There are two reasons for expecting modifications to 
the semiclassical spectrum. The first is related to the 
singular spatial structure of the supervelocity field. One 
flux quantum associated with a single vortex means that 
the superconducting order parameter, or, equivalently, 
the wave function of a Cooper pair (charge 2e) is single 
valued and has a phase winding of 2n around each vor- 
tex. As emphasized in Rcf. the semiclassical approach 
transfers this phase winding equally to each of the quasi- 
particles forming the Cooper pair (each having charge e). 
Consequently, their wave functions change phase by 7r 
around a vortex line ( Aharonov-Bohm phase) , leading to 
the necessity of introducing branch cuts and the problem 
of multi-valued wave functions in the full quantum me- 
chanical treatment. However, the semiclassical approx- 
imation is only valid for large quantum numbers, that 
is for the quasiparticles for which the total phase of the 
wave function, accumulated as the electron moves around 
the vortex, is large. The wave function of an electron 
circling a vortex at a distance r from the vortex center 
acquires a phase 27rfc/r, compared to an extra Aharonov- 
Bohm phase 7r from the supervelocity field. For the anal- 
ysis of the extended states (r 3> £o) m the semiclassical 
approach (valid at fc/£o 3> 1), k/r S> 1, so that if the 
phase of the wave function is changed by n, it still cor- 
responds to the quasiclassical state with essentially the 
same energy and momentum. Since £o ~ w//Ao, we can 
rewrite the condition for the applicability of the semi- 
classical method &/£o > 1 as ao = u//fA 3> 1. In- 
deed, the work of Refs. [35 36 has shown that for large 
anisotropy of the Dirac cone the semiclassical approach 
remains valid down to the lowest energies. Since a ~ 14 
for YBCOeJ, and a ~ 20 for Bi-2212BEl this is the pa- 
rameter range relevant for the study of HTSCs. In a very 
recent preprint Mel'nikov has shown that the Aharonov- 
Bohm phase leads to a different result for the quasipar- 
ticle density at distances r » Ai, whila-in the range 
Co *C r <C Xl the semiclassical results hold.E3 Once again, 
since in the field range where most experimental measure- 
ments are done the intervortex distance d <C Xl this re- 
sult suggests that the semiclassical approach is adequate 
for the analysis of these experiments. 

Quantum mechanical treatment is nevertheless needed 
for accurate descrip tion pf the states at very low energies. 
Kopnin and VoloviktaO have considered the effect of the 
magnetic field on the nodal quasiparticles perturbatively, 
and found that the spacing between quantum mechani- 
cal levels of the near nodal quasiparticles for which the 
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spatial extent of the wave function is comparable to the 
intervortex distance is Ekv = v^/d = Eh I old- There- 
fore they have argued that below this energy scale the 
semiclassical approach becomes invalid. For a ~ 15 this 
energy scale is of the order of a few Kelvin per square 
root of Tesla. However the specific heat measurements 
show no crossover to a novel behavior at that scaleQ, and 
the measurements of the thermal conductivity below 0.5 
K in fields of up to 8 Tesla are in agreement with the 
semiclassical calculations.! 

Marinelli et alSB have suggested on the basis of their 
numerical work that the actual crossover energy de- 
creases much more rapidly with the increasing anisotropy 
a than the result of Ref. |39| implies. They have noted that 
the dispersion of the quasiparticles along the Fermi sur- 
face is strongly softened in the vortex state, so that the 
effective anisotropyL£l an increases much more rapidly 
than ctD, leading to a much smaller crossover scale, per- 
haps exponentially small in a^.cHI Since in real samples 
the presence of impurity scattering and the disorder in 
the vortex lattice always smear out the energy struc- 
ture on small scales, we therefore expect that for the 
purposes of comparison with the measurements of the 
thermodynamic quantities, the semiclassical description 
is adequate. 

Therefore for the parameter range relevant to the study 
of most real unconventional superconductors, the semi- 
classical approach reproduces the energy spectrum of the 
near-nodal quasiparticles in the vortex state to a high de- 
gree of accuracy. Moreover, presently it remains the only 
approach which is capable of including the effect of impu- 
rity scattering into the analysis, and we use it hereafter. 



III. SEMICLASSICAL APPROACH TO THE 
VORTEX STATE 

A. Nodal approximation 

The semiclassical approximation takes as its starting 
point a Fermi-liquid description of the nodal quasiparti- 
cles, so that in the absence of a magnetic field the Green's 
function in the particlc-holc (Nambu) space is given by 
the Bardeen-Cooper-Schrieffer form with an anisotropic 
gap, 



G(k, cjj, 



im n T Q + A k Tj + ( k T 3 

+ G + ■ 



(1) 



Here % for i — . . . 3 are the Pauli matrices (to is the unit 
matrix), ui n — irT(2n + 1) is the Matsubara frequency, 
and £k is the energy of a quasiparticle with momentum k 
measured relative to the chemical potential. We consider 
a two-dimensional Fermi surface with an energy gap of 
d x 2_ y 2 symmetry given by At = Ag(k 2 — k 2 )/k 2 . Low 
energy properties depend only on the nodal quasiparti- 
cles, and are only functions of the parameters entering 
the linearized dispersion near nodes at position k„. As 



Ck ~ v / ■ (k— k„), and Ak ~ va ■ (k— k n ) near a node, the 
poles of the Green's function after analytic continuation 
to the real axis, iui n — > ui + id, are located at energies 



E(k) = ±^a+Ai*± 



jkl 



,,2 1.2 

J A K \\ ' 



(2) 



where k± and ku are the components of k — k„ normal to 
and along the Fermi surface respectively. We parameter- 
ize the Fermi surface near each of the four nodes not by 
the momenta k± and fey , but by the quasiparticle energy 
E and the angle 8 defined as 



Vfkx = E sin 9, 
VAk\\ = EcosQ. 



(3) 
(4) 



The energy cutoff is chosen to preserve the volume of the 
Brillouin Zone, so that for a square latticp-svith the pe- 
riodicity a, it is set at Eq = ^/mJJWK/ aE^ By making 
this choice we extend the conical dispersion law beyond 
the maximal gap value, Ao- This leads to logarithmic, in 
Eq/Aq, corrections to the quantities which depend on the 
cutoff energy. Since i>a ~ A /fcy and kf ~ 7r/a, we ob- 
tain Eo ~ WEfAo, where Ef is the Fermi energy. In the 
high-T c materials Ef ~ 3— 10Ao, and, consequently, the 
choice of Eq as the cutoff energy does not affect the re- 
sults significantly. Therefore near each node the Green's 
function at real frequencies can be written as 



G{E,0;w) 



ojtq + E cos Qt\ + E sin Qt^ 



10' 



E 2 



(5) 



In writing Eq.(|^) we have replaced the bare frequency 
lu by the renormalized frequency uj to include the effect 
of impurity scattering. We account for isotropic strong 
(phase shift tt/2) impurity scattering in the framework of 
a self-consistent T-matrix approximation, and consider a 
particle-hole symmetric system, so that the only non- 
vanisbing component of the self energy is proportional 
to To.c3 Therefore the effect of impurities is to replace in 
the Green's function to by its renormalized value, uj = 
uj — with the self-consistency condition 



E(25) 



(6) 



where u, is the impurity concentration. In the nodal 
approximation the integral over the Brillouin Zone can 
be written as a sum over the nodal regions 
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Writing lu = u>\ + iu)i, we obtain for Eq 3> M 



2 u>\ + iu)2 

7T VfVA 
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(8) 
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The well-known relationships for the density of states in 
the pure limit (w 2 = 0, u>\ — u>), and for the residual den- 
sity of states in the presence of impurities {uj\ = 0, 0J2 = 
7) follow easily (cf. Ref. [ll]) from 



N(w) 



-- V ]lmGn(k,uj), 

TT • ■ 



to give 



N(uo) 
N(0) 



TTVfV A 

2 7 



1 E ° 
2 ln — ' 

TT VfV& 7 



The self-consistency condition oj\ 
the latter case is (cf. Ref. fl3]) 



pure limit 
unitarity. 

iuo 2 



(9) 

(10) 
(11) 

u> — S(5) for 
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B. Doppler shift 

In the semiclassical approach to the the vortex state 
the presence of a superfiow is accounted for by introdue.- 
ing the Doppler shift into the energy uj — > oj + e(k, r)EJ'll3, 
where 



e(k,r) = v s (r) • k, 



(13) 



and v s (r) is the supervelocity field at a position r due 
to all vortices. It was demonstrated by Kiibert and 
HirschfcldEJ that to very high accuracy the Doppler shift 
at the node k„ can be used to approximate the Doppler 
shift for the entire nodal region. Therefore the Green's 
function near each node can be written as 



g(E,@;uj;T)=G(E,@;LO + e n (r)), 



(14) 



where G is given by Eq.(||), n labels the nodes, and 
e n (r) = Vs (r) • k„ . For a d-wave superconductor there are 
two pairs of nodes such that ki = — k 3 and k 2 = — k 4 , 
so that the possible values for the Doppler shift are ±ei 
and ±62- 

In principle now all the physical quantities can be com- 
puted with the help of this Green's function. Several 
comments have to be made about the assumptions im- 
plicitly present in such calculations. First, we neglect all 
inelastic processes. Second, there is no additional quasi- 
particle damping due to the presence of the vortices: in 
the absence of impurities the lifetime of a quasiparti- 
cle defined by the pole of the Green's function in Eq. 
( |l4| ) is infinite. The scattering of the nodal quasiparti- 
cles by vortices depends strongly on the nature of the 
vortex cores; in-|tbc high-T c materials this is an unre- 
solved problemerea If the point of view is taken that 



the core is identical to that of a BCS-like superconduc- 
tor, the neglect of vortex scattering is reasonable for a 
vortex lattice with long range order, and it remains valid 
when only short range order exists provided that the life- 
time is restricted by the impurity scattering rather than 
vortex disorder. The role of disorder in the vortex lattice 
is especially important for transport properties, where 
there is a competition between the increase in the den- 
sity of quasiparticles and the change in the transport 
lifetime in an applied, feld; , several authors have ana- 
lyzed its ccw 
unresolvedc 



•in v>. ael@ and some questions remain 
c3. STM measurements suggest that the 
vortex lattice is ordered in YBCOc3, and that short range 
order is present in Bi-22120, therefore in that case the 
assumption is justified. We note that thermodynamic 
quantities, such as the density of states (DOS), depend 
on a single energy scale, Eh = Vf/d even in a disordered 
vortex state in absence of strong pinning, and this depen- 
dence appears to be nearly-identical for the ordered and 
disordered vortex lattices^. Therefore we expect that 
the results obtained within the semiclassical approach 
remain at least qualitatively valid even for a strongly 
disordered lattice. 

Third, in the analysis of the impurity scattering this 
approach assumes that the positions of vortices and of 
impurities are uncorrelated. The self energy given by 
Eq.(|J) is obtained after averaging over the positions of 
impurities, and solving this equation with the Doppler 
shift included in the Green's function implies that the 
impurity average and the spatial average are taken inde- 
pendently. 

Finally, in the discussion so far we have neglected the 
Zeeman splitting altogether; this is justified when the 
Doppler energy scale exceeds the Zeeman shift. In the 
absence of spin-orbit coupling the Zeeman shift is jiH sa 
0.677? K- T _1 , while the Doppler shift is E H ~ 30VH 
K- T -1 / 2 ; consequently the two become comparable only 
at H cross r~j 10 3 T, and hence the Zeeman splitting is 
irrelevant. On the other hand, for the field applied in 
the plana, the coefficient in the Doppler shift is much 
reducedEfj, and the Zeemaa splitting is relevant for some 
experimental geometries. tB Consequently, we will revisit 
this question in the analysis for this configuration. 

If we know how to express a physical quantity F in 
terms of the Green's function we can now compute its 
local value F(r) with the local Green's function given 
by Eq. (|l4|) . We then approximate the field-deper 
measured value F(H) by the spatial average of F(r)l 1,: 



F ( H ) = J J rf 2 rF(ei(r),e 2 (r)), 



(15) 



where the integral is taken over the part of a unit cell of 
the vortex lattice (with the area A) in real space where 
the Doppler shift is much smaller than the gap maximum. 
Therefore the integration is to be cut off at distances of 
the order of £0 from the center of each vortex. In prac- 
tice in many cases the contribution of the core region 



5 



( r < £o) is small due to the geometric effect (integrals 
are weighted with the surface area rdr) and the integral 
can be extended to the entire unit cell. We note that the 
averaging procedure is often non-trivial for response func- 
tions; for the thermal conductivity K, for example, «(r) or 
1/k(t) are averaged depending on the relativ^pwentation 
of the magnetic field and the heat current .E-jED The av- 
erage in Eq.(p^|) depends on the distribution of vortices. 
In practice, this spatial average has been computed an- 
alytically only for the supervelocity field corresponding 
to an isolate*!-. fjk*x line, cut off at the average intervor- 
tex distance ,t3~t3 and numerically for the pancake liquid 
state.ca 

The starting point of our approach, which simplifies 
calculations and makes possible a generalization of the 
semiclassical method to an arbitrary configuration of vor- 
tices is to rewrite the average as the integral over the 
probability distribution of the Doppler shift for a particu- 
lar vortex configuration. There are, in general, two types 
of local quantities, and therefore of averaging procedures, 
which are required. The density of states in the absence 
of impurity scattering, for example, is a direct sum of the 
contributions from each node, 



AT(o; ) r) = -i-Im|^Trg(k ) a;)| 



(16) 



and can consequently be expressed as an integral over the 
probability density of the Doppler shift at a single node, 

N(lu,H) = ±J2 deN(u; + ae)V(e), (17) 

1 a = ± J-oc 



where 



V(e) = J J ^r^e-v 5 .(r)-k, 



(18) 



Such an approach has been recently used to analyze the 
behavior r©f the interlayer conductivity in the vortex liq- 
uid statecS, where the function V was determined from 
numerical simulations. A similar method (although with 
an unrealistic distribution, see below) has, been used in 
the analysis of the thermal conductivityOo. 

However, in general the function F depends on the 
values for the Doppler shift at two inequivalent nodes, t\ 
and e 2 , and the corresponding average can be written as 

/+oo 
de 1 de 2 F(e 1 ,e 2 )C(e 1 ,e 2 ), (19) 
-OO 

£(ei,e 2 ) = i y"d 2 r^ £l -v s (r) ■ (20) 
x6[ e 2 - v s (r) ■ k 2 ), 



where ki and k2 label two nearest nodes. This is the 
case, for example, for the density of states in the presence 
of impurity scattering, since the self energy (implicitly 
present in the Green's function in Eq. (|l6|)) contains the 
sum over the nodes, see Eq.(^), and therefore depends on 
both ei and e 2 . In general, the function C has to be even 
in both ei and e 2 , and symmetric under the interchange 
(\ £2! in all the cases considered below it depends on 
a single variable t\ + e|. 

Now all the relevant information about the structure of 
the vortex state is contained in the functions C(e\, e 2 ) — 
C'{e\ + el) and 



T(e) 



J rfei£(e, ei), 



(21) 



and therefore to analyze the field dependence of the phys- 
ical quantities we first focus on determining these prob- 
ability densities. 



IV. PROBABILITY DENSITY FOR THE 
DOPPLER SHIFT 

The distributions V and C can be determined numer- 
ically for an arbitrary configuration of vortices. Here we 
are interested in making progress analytically, and there- 
fore consider several model configurations for which the 
distributions can be found exactly. Moreover, we propose 
that the distributions which we consider give the maxi- 
mal and the minimal possible weight to the low-energy 
Doppler shift, and therefore can be used to obtain the up- 
per and the lower limits of the experimentally accessible 
quantities. 



A. Single vortex, H||c 

The simplest of these models is that of a velocity field 
of an isolated vortex, cut off at the distance equal to 
the intervortex distance; since the experiments are in 
the dilute vortex limit such an approach gives an ade- 
quate description of the vortex state. The supervelocity 
is v s (r) = Ti6/2mr, where 9 is the winding angle of the 
vortex in real space, m is the effective mass, and r is the 
distance from the center of the vortex. We now write the 
Doppler shift in terms of the energy scale Eh = Vf/(2R), 
where R is the radius of the unit cell 1 of r the vortex lattice, 



taken to be circular, R = $0/71- I 



v s (r) • k f 



hkj- 
2mr 



sin# 



E 



H 



sin 9. 



(22) 



Here we have introduced the normalized length p = r/R, 
and have chosen, without loss of generality, k„ along the 
direction 9 = ir/2. 

The probability distribution at a single node is now 
easily obtained from Eq.(|l8|) 
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V(e) = - [ d6 j pdps(e-—sint 



(23) L(x,y) = Ej I C(e 1 ,e 2 ) = 



1 



El 



dO *nr (){-)(' ' 



^ sinflje^l- — sin(9 j. 



yielding 



1 E- 
7r e' 



H 



arcsm ■ 



if e> E H ; 
if e < E H - 

(24) 

Here we have taken e > 0, the probability density is even 
in e. 

It was argued in Ref. ^ that the function V{e) for 
any vortex configuration has two important properties. 
First, the asymptotic behavior V(e) = E H /(2e 3 ) holds 
for A » e ^> Since the vortices repel each other, 
the vortex cores do not overlap. The large Doppler shifts 
come from the regions near the cores, where the super- 
fluid velocity is high, and consequently are dominated 
by the single vortex physics. Second, in the absence of 
strong pinning V{e) has a single energy scale E^ and 
depends on the Doppler shift only via c/Eh- Since the 
probability density is normalized, 



/+oo 
V(e)de = 1, 
-OO 



(25) 



we can follow Ref. and define a normalized dimension- 
less probability density as 

P(x) = E H V(e/E H ), (26) 

where x = e/En- 

The two-node probability distribution function is 



C(e u e 2 ) = - [ del 
k Jo Jo 



/*//**(<:!-— sin 0) (27) 



X S I e 2 - 



d sin(6> + 4> ) 
sva.9 



where </>o is the angle between the nodes ki and k2 at 
the Fermi surface. For the pure d-w&ve symmetry which 
we consider here, 0o = tt/2, and the integral can be eval- 
uated to give 

£(ei,e 2 ) = ~ ? 2 E K, 2 , XV4 + 4>E H , (28) 

and zero otherwise. The physical reason for the discon- 
tinuity is that for the nodes at the orthogonal positions 
( \ + e 2 = E\l p 2 > E^, so that the probability of hav- 
ing the Doppler shifts not satisfying this inequality is 
identically zero. In an orthorhombic system, where the 
nodes are not at angle tt/2, the shape of the distribution 
is different. In analogy with the single node probabil- 
ity density we can also define the dimensionless energies 
(x,y) = (ei,e2)/En, and introduce the function 



7r(x 2 + y 2 )" 2 



if x 2 + y 2 > 1, 
(29) 



and zero otherwise. 



B. Single vortex, H| ab 

For the magnetic field applied in the superconducting 
plane it has been recently argued that for a relatively 3- 
dimensional high-T c material, such as YBCO, the semi- 
classical approach still captures the essential features of 
the quasiparticle behavior Jia The approach of Ref. [l6| is 
to take the supervelocity field from an anisotropic Lon- 
don model, but to introduce the Doppler shift only in the 
dispersion of the quasiparticles with the momenta in the 
plane. After rescaling the c-axis to make the unit cell 
of ±he vortex lattice isotropic, the Doppler shift is given 
byH 



v s (r) • kf — sin#sin(0 — a), 



(30) 



where the angle (f> parameterizes the cylindrical Fermi 
surface, a is the angle between the direction of the mag- 
netic field in the plane and the x-axis, and the in-plane 
energy scale is E a b = rjEu, where in the London effec- 
tive mass model the anisotropy rj = {\ a b/ ^c) 1 ^ 2 ■ In the 
nodal approximation (which proiddes an excellent agree- 
ment with the numerical resultstj) the probability distri- 
bution of the Doppler shift at a single node is given by 
Eq.(f24j) with E H replaced by E\ = E ab \ sin(7r/4 — a)| and 
E 2 = E a b \ cos(7r/4 — a) | respectively for the two pairs of 
nodes. Any effects of the three diaiensionality reduce the 
effective value E a b rather severelyc3, so that the estimate 
obtained using the value of rj for the effective anisotropy 
in the two-dimensional case can only serve as an upper 
limit. 

For such a geometry the Doppler shifts at the two 
neighboring nodes are related by e 2 — E 2 ei/Ei; in con- 
trast to the case of the field applied along the c-axis, the 
Doppler shift at one of the nodes uniquely determines 
the value of the Doppler shift at the other node indepen- 
dently of the winding angle 9 in real space. Therefore 
the two-node probability distribution is given by 



(31) 



and a single average is always sufficient for computing the 
physical quantities in the semiclassical approximations 
for the field applied in the plane. 



C. Vortex solids and liquids 

We now discuss how the probability densities obtained 
above can be generalized to the case of vortex solids or 
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liquids. We first consider the single node probability den- 
sity P{x). Since this function is normalized, the question 
is what type of the redistribution of the density in Fig.l 
one may expect for realistic vortex structures. As argued 
above, the high energy tail of the distribution is entirely 
determined by the single vortex physics, and is therefore 
insensitive to the structure of the vortex state; the redis- 
tribution of weight occurs in the region x < 1 or e < Eh- 

It is also clear that the single vortex picture described 
above underestimates the number of points where the 
Doppler shift vanishes. For the supervelocity field of a 
single vortex |v s (r)| > everywhere in the unit cell, and 
the Doppler shift vanishes only for the superfluid velocity 
direction normal to the nodal directions in k-space. In 
a vortex lattice there exist points where |v s (r)| = 0: the 
high symmetry locations such as midpoints between the 
centers of two neighboring vortices. Consequently, for 
vortex lattices -P(O) is larger than it is in the single vortex 
picture. The weight shifted to the vanishing Doppler 
shift, x = 0, comes at the price of a reduction in the peak 
in P(x) and moving the peak to smaller x. The actual 
shape of the function depends on the type of the vortex 
lattice, the number of nearest neighbors, and on relative 
orientation of the basis vectors of the vortex lattice with 
respect to the nodal directions. 

As the number of nearest neighbors is increased, so is 
the value P(0). This value depends not only on the num- 
ber of zeroes, but also on the asymptotic behavior of the 
supervelocity near a point where it vanishes, v s (tq) = 0. 
It is easy to check that if v s (r — r ) oc |r — r | T ', the contri- 
bution of this area to P(0) is finite for 77 < 2, is singular 
but integrable for 2 < rj < 3, and is non-integrable (and 
therefore non-physical) for 77 > 3. In a typical vortex 
distribution v s varies linearly with the distance from ro, 
so that P(0) remains finite. We now try to derive ana- 
lytically an approximate distribution which gives a large 
weight to the probability of the vanishing Doppler shift; 
we consider it here to model a relatively disordered vor- 
tex state, such as a vortex liquid, and to provide a lower 
limit of the magnetic field dependence of the physical 
quantities. To make progress we consider a cylindrically 
symmetric spatial dependence of the supervelocity, mod- 
ulated compared to the single vortex distribution. Differ- 
ent choices for the modulation-of the superfluid velocity 
are considered in the literatureEIrE3; in any approach the 
supervelocity near the vortex core should remain nearly 
unmodified compared to the single vortex velocity field, 
while at the cell boundary v s = 0. Therefore, in the 
cylindrically symmetric case, the Doppler shift (for H||c) 
can be approximated as 



v s (r)-lL f = E H S(p)sm6, 



(32) 



vortex distribution), the required asymptotic behavior of 
S(p) is different from that of v s in the system with points 
of vanishing Doppler shift. Nevertheless, as we show be- 
low, the appropriate choice of S(p) allows us to arrive at 
a probability distribution close to that obtained by nu- 
merical simulations of the vortex liquid. In such a liquid 
the distribution is temperature dependent. A detailed 
calculation therefore would have to take into account the 
changes in the probability density with the temperature 
in a given material. These changes are not well under- 
stood beyond simple models, and even then are usually 
accessible only via a dynamical simulations. We there- 
fore take the point of view that for a qualitative or semi- 
quantitative analysis it is sufficient tOj-cpnsider a model 
temperature-independent distribution.^ 

Computing the distribution V(e) from Eq.(|lS|) we ob- 
tain 



P{x) = - [ 
* Jo 



pdp 



VsHpT 



(33) 



Clearly, P(0) is finite when S(p — ► 1) oc (1 — p) n with 
Tj < 1. We use here two different models where the 
superfluid velocity field of a single vortex is modulated 
to vanish at the unit cell boundary in a fashion which 
allows analytical progress. In the first, we take the 
modulating factor to be (1 — r 2 / R 2 ) 1 / 2 , which leads to 
S{p) = \/l — p 2 1 P- Computing the probability densities 
as in the previous section we find 



£(ei,e 2 ) 



1 



El 



n(el + 4 + E%r 



and 



P(e) = 77 



Eh 



2 (E 2 H + e 2 ) 3 / 2 ■ 



(34) 



(35) 



In this case the measure of Doppler shift zeroes is large 
due to disorder in the positions of vortices, P(0) = 0.5. 
The simplicity of this probability distribution makes this 
choice attractive for further analytical work. 

Another possible choice is S(p) = y/1 — p/p; it leads 

to 



£(ei,£2) 



1 



Eh 



c ?\3 



(36) 



E 2 

H 



1 



E 2 

H 



Note that as (e 2 + t%)Eh ~ * tne distribution L is fi- 
nite: £(0,0) = 2/(wE H ). The corresponding single node 
probability density is given by 



where S(p -► 0) oc l/p and 5(1) = 0. 

Notice that the requirement that P(0) is finite imposes 
restrictions on the decay of S(p) as p — > 0. Since in the 
cylindrically symmetric model v s vanishes along a line 
rather than at discrete points (as it does for a realistic 



P(e) 



1 



irE 



H 



Eh 



3Eh 
4e 5 



1 



arccos 
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For this distribution P(0) = 32/[157r) « 0.68, larger than 
the value of 0.5 given by Eq. (|35|). 

The probability density P(x) for all three distributions 
is shown in Fig. 1. In the following we will refer to the 
distributions given by Eqs. (j34)-{ 35 ) and by Eqs. (p6|)- 
([57]) as liquid I and liquid II respectively. The reason for 
that is clear from the inset of Fig.l: these distributions 
are close to those obtained with the help of the Langevin 



dynamics simulations of the pancake liquid in Ref. 48 



in the vortex liquid they preserve the cylindrical symme- 
try of the supervelocity field on average, while introduc- 
ing zeroes in that field because of the cancellation of the 
supervelocity from neighboring vortices. For a realistic 
vortex lattice we expect the results for thermodynamic 
quantities to be bracketed by the values obtained in the 
single vortex approach, which overestimated the effect 
of the field by under-counting the number of points in 
the unit cell of the vortex lattice where the Doppler shift 
for quasiparticles near a particular node vanishes, and, at 
least approximately, by the liquid II distribution given by 
Eqs. (p6|)-(|3~7|). The distribution function for the pancake 
liquid can be even sharper peaked at x = 0, nevertheless 
we believe that the approximate analytic form provides 
a reasonable low-end estimate for most experimental sit- 
uations. 



Q_ 




FIG. 1. Main panel: probability distribution P(x) in the 
single vortex approximation from Eq. ( |24[ ) (s olid line), and 
for a model vortex liquid states from Eq. feq) (dashed line) 
for liquid I model, and from Eq. ( pS7| ) for the liquid II model 
(dot-dashed line). Inset: comparison of the distributions for 
the model liquid states (same notations as in the main panel) 
with the numerically determined distributions for pancake liq- 
uid in BSCCO at T = 9K (narrow distribution) and T = 67K 
(broad distribution) from Ref.48. 



D. Completely disordered vortex state 

The "universal" high Doppler shift behavior of the 
probability density, P(x) cx a; -3 results from the strong 
repulsion between the vortex lines which prevents vor- 
tex cores from overlapping. If the vortices were non- 
interacting, in a disordered state their positions would 
be completely random, leading to a Gaussian distribu- 
tion of the Doppler shifts. Such an .-approximation has 
been used by Yu et aZx3 and FranzeS in their analysis 
of the thermal conductivity in the vortex state. Even 
though it is never realized, it is instructive to compare 
the predictions obtained with such a distribution with the 
results obtained in the framework outlined above. The 
comparison may be useful for the extremely anisotropic 
layered superconductors in the geometry with the field 
applied in the basal plane. In that arrangement vortices 
lack proper cores, the intervortex repulsion is weakened, 
and we expect significant disorder in vortex positions due 
to the presence of defects (such as boundary effects, twin 
boundaries, etc. ). Consequently, the 1/x 3 asymptotic 
behavior does not onset up to large Doppler shifts (very 
close to the core), and over the low (compared to the 
gap amplitude) energies, the probability density decays 
rapidly. We therefore also consider in the following the 
random distribution of vortices, which leads (omitting 
factors of In /£q ria-the width of the Gaussian) to the 
probability densityEylla 



P(x) 



1 



(38) 



We now investigate the dependence of the thermody- 
namic coefficients on the magnetic field and the temper- 
ature for different structure of the vortex state and com- 
pare it with the experimentally observed behavior. 



V. DENSITY OF STATES: PURE LIMIT 

We begin by considering the density of states and the 
electronic contribution to the specific heat in the pure 
limit. While this is one of the simplest quantities to ana- 
lyze, it is the one directly relevant to the measurements of 
the fieldjdependence of the specific heat in YBCO single 
crystals.! 2 ]" u To justify ignoring impurities in this analysis 
we emphasize that the energy scales associated with the 
Doppler shift are quite large, and at moderate fields ex- 
ceed the impurity bandwidth even in not too clean ss 
pies, and exceed it by far in the latest single cryst 
Taking the Fermi velocity Vf ~ 1.5 — 2.5 x 10 7 cm/sE2l, we 
obtain E H /V~H ~ 30 K-T" 1 / 2 , and for YBCO near opti- 
mal doping, where I/77 ~ 2.5—4, we obtain E a f,/VH < 10 
K-T- 1 / 2 , while the impurity bandwidth 7 is of the order 
of a few Kelvin or less. 







A. Density of states for H||c. 

We first consider the experimental arrangement with 
H||c. The density of states in the pure limit is given by 
Eq.@ leading to 



N(w,H) 



deN(u + e)V{e) 

TTVfVA 



(39) 



Introducing the dimensionless variable x — e/Eu and 
considering hereafter u> > we find that the density of 
states is given by 



2b) i-u/Eh 
N(u,H) = / P(x)dx 



irvfVA Jo 



(40) 



2E 



H 



TTVfVA J^/Eh 



xP(x)dx. 



The scaling properties of the density of states with 
u>/EjjE2 can be made obvious by rewriting it as 



TTVfVA ^Eh' 



(41) 



F N {Z) = 2^Z P{x)dx + J xP{x)dx^j. (42) 

The residual density of states at the Fermi surface is 
given by 



TTVfVA <Zva V ^fPo 



H 



(43) 



where M\ is the first moment of the probability distribu- 
tion of the Dopplcr shift 



Mi = 2 / xP(x)dx, 
la 



(44) 



which contains all the information about the structure 
of the vortex state relevant to the magnitude of the 
\[H term in the specific heat. For the probability den- 
sity given by Eq.(p3) (single vortex model) we then find 
Mf ~ A/tt k, 1.27, while for the liquid I distribution 
given by Eq.(^) we obtain M\ = 1. For liquid II distri- 
bution the integral can be evaluated numerically to give 
M 12 ps 0.85, while for the completely disordered distribu- 
tion of vortices Mf = l/y/n ps 0.56. We therefore expect 
that M\ ~ 1 for any realistic vortex state. Furthermore, 
since the number of zeros of the Dopplex,shift increases 
with the increased disorder in the latticeo, we expect on 
general grounds that the coefficient is larger for the more 
ordered vortex state. The residual density of states given 
by Eq. ( |43"| ) is close to the expression obtained, by Won 
and Maki in a different approximation schemd£3. 

Expanding Eqs.(|4l])-(p2|) at low energies lo <C Eh we 
find 



N(u>,H) 



E 



H 



VfVA 



Mi 



P(0) 



(45) 



Therefore the energy dependence of the density of states 
in the field dominated regime is determined by the prob- 
ability weight of the vanishing Doppler shift. As the lat- 
tice changes toward a larger coordination number and 
towards disorder, the measure of points where the su- 
perfluid velocity vanishes increases. As a result, the co- 
efficient of the leading field dependent term, Mi, de- 
creases, while the coefficient of the energy dependent 
term, P(0), increases. The values of this coefficient are 
P s (0) = 2/3vr re 0.21, P l (0) = 0.5, P l2 {0) = 0.68, and 
P 9 (0) — l/v^ ~ 0-56 for the single vortex, liquid, and 
Gaussian distributions respectively. It immediately fol- 
lows that the position of the crossover from the field dom- 
inated to the zero-field temperature dominated behavior 
in the average density of states is much more sensitive to 
the structure of the vortex state than the leading field- 
dependent term. 

In the effective weak- field range, u> ^> Eh, the field de- 
pendent contribution is independent of the distribution 
of vortices. The vortices are well separated, and the re- 
gions where the Doppler shift exceeds the temperature 
are close to the cores, and consequently dominated by 
the universal tails P(x) = l/2x 3 , yielding 



N(w,H) 



TTVfVA 



IE 



1 + -- 



H 



2 cj 2 



(46) 



The full dependence of the density of states on the 
ener gy and the magnetic field can be obtained from 
Eqs.(^l|)-(^2|) with the probability densities discussed 
above. For the single vortex_picture we regain the re- 
sult of Kubert and HirschfeldEj 



F* N (Z) = - 

TT 



TT l + Z-72 



(47) 
if Z > 1; 



(1 + 2Z 2 ) arcsinZ + 3Z^1 - Z 2 , if Z < 1 



For the liquid I model we obtain a remarkably simple 
result 



F l N = y/Z 2 + 1, 



N 1 (lo,H) 



Vu; 2 + E 2 H 

TTVfVA 



(48) 
(49) 



while for the liquid II model the integral can only be 
evaluated numerically, and for the Gaussian model 



F%{Z) = Z^{Z) 



exp(-Z 2 ) 



(50) 



where $>(Z) is the probability integral.o Notice that for 
the Gaussian distribution the enhancement of the density 
of states in the weak field limit, to ^> Eh is vanishingly 
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small in uj/Eh, in contrast to Eg. (fl6|) . Indeed, in this 
limit the field-dependent part of the density of states is 
determined by the weight in the part of the distribution 
V(e) with e > u>, which is exponentially small. 




1.0 2.0 

Q)/E H 

FIG. 2. Energy dependence of the density of states in the 
magnetic field for different models of the probability den- 
sity for the Doppler shift. Density of states is in units of 
-Eh/(tto/«a). 

This difference is clear from Fig. 2. The low energy 
limit of the density of states depends on the moment of 
the distribution function, and is therefore different for 
each of the model distributions. On the other hand the 
high energy, or weak field, limit yields the same result for 
the models respecting the asymptotic a; -3 decay for the 
probability distribution P{x), while the Gaussian model 
gives the density of states which is not enhanced relative 
to the zero-field value. The Gaussian model therefore 
misses the field-dependent contribution to the physical 
quantities at high energies, leading to incorrect results, 
especially in the regime T < Eh- 



B. Density of states for H| ab. 

We can now analyze in the same framework the 
anisotropy in the density of states for the field applied 
in the superconducting plane . at an angle a to the x- 
axis. As discussed in Section [VB , density of states for 
such a configuration is a sum over the two inequivalent 
pairs of nodes, with the different characteristic scales for 
the Doppler shift at each pair of nodes, 



Ei = £ ah |sin(7r/4- a)|, 
E 2 = E o6 |cos(7r/4-a)|. 



(51) 
(52) 



In the London model E ab — r/En, where r\ is the pen- 
etration depth anisotropy ratio. We emphasize that in 
reality the value of the "effective" anisotropy depends 



on the details of the c-axis transport propertiesB, and 
therefore the estimate of E ab /\f~H ~ 10 K-T -1 / 2 is just 
an upper limit on its magnitude, and, as we comment 
below, the value inferred from the available experimental 
data on the specific heat is lower. 

The density of states in the clean limit is given by 



N{uj,H-q) 



1 



N 1 (lo,H) + N 2 {uj,H) 



(53) 



where N% is computed from Eq.(^Tj) as in the previous 
section but wifckp& (i = 1,2) replacing Eh- 

As is knowrilaH the residual density of states exhibits 
fourfold oscillations as a function of the direction of the 
applied field in the plane 



N(0,H; a) 



Mi E l + E 2 

2 TTVfVA 

Mi E ab 



(54) 



V27T VfVA 



max since , cosa 



The minima of the density of states occur when the field 
is along the nodal direction, a — n / A + nn/2. In that case 
at two of the four nodes the circulating currents are in the 
plane orthogonal to the direction of the Fermi momentum 
at the node, and consequently the Doppler shift vanishes 
at all points in real space (either E\ = or E 2 = 0), 
as seen in Fig. 3. In contrast, when the field is along the 
antinodal direction, the Doppler shift in non-zero, and 
all four nodes contribute tcj-the density of states, leading 
to a maximum in N(u, H)t3 




2/ 




H 




3 

\ 



FIG. 3. Contribution of different nodes to the density of 
states. The nodes are numbered and the direction of the 
Fermi momentum is shown at each nodal point. Left: field 
along the nodal direction and orthogonal to the other pair 
of nodes. The Fermi momentum at nodes 2 and 4 (broken 
line) is orthogonal to the plane where supercurrents flow, and 
the Doppler shift vanishes everywhere in space for this pair of 
nodes. Right: field in the antinodal direction. The Doppler 
shift is non-vanishing at some points in space for each node, 
and the density of states is maximal. 

It is important to emphasize that, as is clear from 
Fig. 3, it is only for a tetragonal system that the min- 
ima in the density of states occur for the field along the 
node. One reason for that is that in an orthorhombic sys- 
tem (m a 7^ To;,), for the field applied in a direction other 
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than along the principal axes of the effective mass ten- 
sor, the. directions of the internal and the external fields 
differ £3 The difference may be quite small in the exper- 
imentally relevant field range; ignoring it, Schachingcr 
and CarbotteE3 argued that the minima occur when the 
field is parallel to the direction of the Fermi velocity at 
the node, which differs from the direction towards the 
node. 

The anisotropy in the density of states given by Eq. ( p34| ) 
is ~ 30% for the purely two-dimensional model consid- 
ered here. Any three-dimensionality reduces this number 
severely: if there is a line of nodes extending along the 
z-axis, for the field applied towards a node in the equato- 
rial plane, the Doppler shift vanishes only for the nodal 
quasiparticles with momenta in the plane. For the quasi- 
particles on the same nodal line but with a component 
of the momentum along the z-axis the Doppler shift is 
finite, consequently the node is still "active" in contribut- 
ing to the density of states. In the simplest estimate4ti 
a 3-dimensional system the effect is reduced to 7-8%li3, 
in some models with a tight binding dispersion along the 
c-axis it may be reduced even further, to about 4%E3, 
making the effect more difficult to detect. 

The anisotrpay is also rapidly washed out with in- 
creased energytS. Since the density of states has a min- 
imum when the field is applied along a node (E\ = 
for example), the corresponding pair of nodes is "inac- 
tive" and insensitive to the field; therefore the density of 
states increases linearly in energy, as in the absence of 
a field. For the field away from the nodal direction the 
density of states increases as a square of the energy, see 
Eq. (|45|), resulting in a rapid suppression of the differ- 
ence between the two geometries. For low energies in the 
limit uj <C E\,E-2, which can only happen if the field is 
not close to a nodal direction (E\, E% ^ 0), we have 



N(uj,H;a) w max[| sin a\ , | cos a\] 
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FIG. 4. Angular dependence of the density of states, mea- 
sured in units of E a i/{-KVfVA), on the the direction of the 
applied magnetic field. Density of states has been computed 
with the model liquid I probability density of the Doppler 
shift. Angle a is measured with respect to the x-axis, and the 
minima are along the position of the nodes. Notice a signif- 
icant reduction in the anisotropy at energies of the order of 

E a h- 

The angular dependence of the density of states van- 
ishes at higher energies, as for uj ^> Ex,Ez with a re- 
alistic distribution respecting the asymptotic behavior 
P(x) tx x~ 3 for x ^> 1 
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(57) 



The exact crossover scale from the strong to the weak 
field regime depends on the particular choice of the prob- 
ability distribution. This is shown in Fig. 5 for the three 
different choices of P(x) considered in this work. 



On the other hand, if E\ <C uj <C E2, which may happen 
when the field is close to one of the nodes, and E\ <C E 2l 
we have 



N{uj,H;a) 



1 



2jTVfVA 



uj 2 E 2 
M 1 E 2 +uj + —P{Q) + -^ 
E 2 2ui 



(56) 



The anisotropy in the density of states as a function of 
the angle for the liquid I model is shown in Fig. 4, at low 
energies there is no qualitative difference between the 
different models, see below. As the energy is increased 
the sharp minima fill up, and the resulting anisotropy 
decreases. 
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FIG. 5. Relative anisotropy in the density of states, 
5N(w) = N(oj,H;a = 0) - N(uj,H;a = tt/4) normalized 
by the maximal Nmax(w) = N(ui,H;0) for different models 
considered in the text. Notice that the relative anisotropy at 
lu — is identical for all models, as is clear from Eq. (pj). In 
the single vortex model the anisotropy vanishes identically at 
u > E a t, see Eq.(^7j). For the Gaussian model the exponen- 
tial asymptotic behavior of the probability distribution leads 
to the inverse anisotropy in the intermediate energy range. 
The two liquid models yield a very similar dependence of the 
anisotropy on the energy. 



C. Zeeman splitting 

Typically the Zeeman shift is small compared to the 
Doppler energy scale, and does not modify significantly 
the density of states. Indeed, for the field along the c- 
axis, the spin-up and spin-down density of states is given 
by 



N ±^H) = -^-F N[ ,. 

ITVfVA \Eh> 



(58) 



where u)± = \u) ± jiH\. Therefore, for example, the cor- 
rections to the residual density of states due to the para- 
magnetic contributions in the regime /J.H <C Eh are of 
the order 



5N(0,H) _ //ig\ 2 P(0) 
N(0,H) ~ \E H ) Mi 



< 1. 



(59) 



For a quasi three dimensional materials, such as YBCO, 
the relevant energy for the field applied in the plane is 
E a i,. Even if this energy scale is only 15% of Eh, a fac- 
tor of 2-3 smaller than estimated, the ratio E a \,/ /iH ~ 
6.7 /VH T 1 / 2 implies a crossover field of 45 T. Conse- 
quently the Zeeman splitting is unimportant compared to 
the Doppler shift in the experiments performed so far in 
YBCO. This is in contrast to more two-dimensional ma- 
terials, such as BSCCO, where the response to a paca 
magnetic field is dominated by the Zeeman splitting!!' 

The situation is different for the geometry with the 
field applied along a node; this was first pointed out in 
Ref. |l9|. In this case the Doppler shift at one pair of the 
nodes vanishes, and, at u> — 0, the only contribution to 
the density of states at these nodes is due to the Zeeman 
splitting. The Zeeman splitting leads to a finite contri- 
bution to the density of states which is linear in [iH , and 
consequently reduces the residual anisotropy SN(0,H). 
This reduction has been investigated in Ref. |lj| 

This is however not the only effect of the paramagnetic 
coupling. Since the density of states at the nodes with the 
vanishing Doppler shift is now dominated by the Zeeman 
splitting at low energies, the anisotropy is not reduced 
as rapidly with the increasing energy. Indeed, if only the 
paramagnetic effect is take into account, the total (per 
particle, i.e. summed over the spins rather than per spin) 
density of states is 



N z (w,H) 



afiH\ 2 max[/ii/, 



TTVfVA 



ITVfVA 



(60) 



and therefore the anisotropy increases with u> up to u) — 
\xH, where it reaches a maximum. 

Let us consider the liquid I model, where the analytic 
expression for the density of states is particularly simple, 
the results are not modified substantially if other models 
are used. The anisotropy in the total (summed over the 
spin directions) density of states between the nodal and 
the anti-nodal directions is given by 



8N{lj,H) = 



1 



2irvfVA 



(61) 
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As Fig. 6 demonstrates, even though the zero-energy 
anisotropy is severely reduced upon inclusion of the Zee- 
man splitting, the anisotropy at moderate energies is 
close to the result obtained without accounting for the 
paramagnetic effect. It is clear that the magnitude of 
the reduction and the crossover energy depend on the 
actual value of E a b, and we need a realistic estimate of 
this value to evaluate the impact of the paramagnetic 
splitting on the experimental results for the field along a 
nodal direction. Such an estimate can be obtained from 
the analysis of the data on the specific heat which we 
discuss in the next section. 
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ab 

FIG. 6. Effect of the Zeeman splitting on the anisotropy 
in the density of states. Here E ab /\fH = a K-T~ 1,/2 . Main 
panel: energy dependence of the anisotropy for the liquid I 
distribution with (solid line) and without (dashed line) ac- 
counting for the Zeeman splitting for a — 5 at H = 10 tesla. 
Inset: anisotropy in the residual density of states, in units of 
a/2nVfVA, as a function of the applied field. Dashed line: no 
Zeeman effect, dot-dashed line: a = 10, solid line:a = 5. 
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VI. SPECIFIC HEAT AND SCALING 

The information about the density of states is exper- 
imentally available primarily via the specific heat mea- 
surements, and we now address this quantity in more 
detail. The preliminary analysis of some of these issues 
within tho, single vortex picture has been carried out by 
us beforeE3. Here we concentrate on the effects of differ- 
ent distributions, and on the measurability of the specific 
heat anisotropy. 



A. Specific heat 

■-electronic contribution to the specific heat is given 



byt 

C(T,H) 



+ 00 



dwN{u,H)[ — I cosh"' — (62) 



T I dxx 2 N(xT, H) cosh" 
Jo 



Making use of Eq. ( ^l| ) we can rewrite the specific heat in 
the form useful for further analysis. For the field along 
the c-axis, Hllc, we have 



C(T,H) = 



2TE H 



nvfVA Jo 



dxx 2 F N (J^- S ) cosh 2 |, 



(63) 



where Fjy is given by Eq. ((42[) . 

As a result we find in the limit Eh T 

TTVfVA v 3 15 E H 

and in the opposite limit, Eh -C T, 
2T 2 



C(T,H) 



(9C(3) 



||m2 



(64) 



(65) 



Two quantities which can be compared with experi- 
ment are the coefficient of the T 2 term in absence of the 
field, 



7s 



k% nV mol 18C(3) 



(66) 



H S TTVfVA 

and the coefficient of the T^/H term at low temperature 
C(T,H) _ tt 3 / 2 k% nV mol Mi 



p = lim 



T 



Va^q' 



(67) 



along the c-axis, and n is the number of CuC>2 layers per 
unit cell. The presence of both terms has been firmly 
established from the analysis_jQf the experimental data 
on the specific heat in YBCOotrH, however, there remains 



disagreement about the values of the coefficients between 
different groups. 

For that material V mo i = 104.6 cm 3 /mol, n = 2, and 
s rj 12A. The coefficient p can in general be determined 
to a higher degree of accuracy, and the values available 
in the literature are p ss XL91 mJ mol" 1 K~ 2 T" 1 / 2 for 
moderately clean samplesfcm, and more recently obtained 
p w 1.3(1 mJ mol -1 K~ 2 T -1 / 2 for the ultra-pure single 
crystalaj. The analysis of these data in the single vor- 
tex picture has been carried out by Wang et alU, and by 
Chiao et alcA. In that picture the 50% difference in the 
coefficient translates into the same relative difference in 
the value for the slope of the gap. In contrast, according 
to the previous section, the more ordered vortex state 
leads to a larger first moment of the distribution, and 
consequently to a larger value of p in Eq.(|67|); it is there- 
fore reasonable that a higher quality crystal would have a 
more ordered vortex state and hence a larger coefficient p. 
If we set Mi = 1 the experimental values of p lead to the 
values for the slope of the gap of va ~ 1-5 x 10 6 cm/s and 
va ~ 1.0 x 10 6 cm/s respectively. On the other hand, 
taking Mi = A /it for the pure crystal, yields a larger 
va ~ 1-27 x 10 6 cm/s, leading to a less than 20% discrep- 
ancy between the groups. The disagreement can be fur- 
ther reduced by assuming a disordered state with Mi < 1 
in the ceramic sample of Ref. We also note that the 
pure crystal of Ref. @ is overdoped, rather than optimally 
doped as in the work of Ref. g,^, which may contribute 
to the difference in the coefficient. In combination with 
the value for the ratio vj/va ~ 14 obtained from the 
universal limit of the thermal conductivity^ this yields 
Vf ~ 1.8 x 10 7 cm/s. This is in reasonable agreement with 
the value of the Fermi velocity obtained from the ARPES 
measurements in BSCCOj 3 !] which is believed to have a 
Fermi surface similar to that of YBCO. 

The coefficient 7 S of the temperature dependence has 
been measured with significantly larger error bars, and 
the results from different groups vary significantly: Moler 
et alu reportecLthe value of 0.1 mJ mol -1 K~ 3 , Wright 
and co-workerso obtained j s ~ 0.064 in the same units, 
while Wang et ala measured 0.21. From the compari- 
son with Eq.(66) we find VfVA ~ a x 10 13 cm 2 /s 2 , where 
a = 2.9,4.5,1.4 for the three values given above. All 
these yield the Fermi velocity within a factor of two of 
the estimate given above. This implies that in the cal- 
culations requiring a cut-off in energy the cut-off Eq ~ 
1.3-2.3 x 10 3 K. 

We now turn our attention to the field applied in the 
a6-plane, and discuss the specific heatJbllowing the gen- 
eral approach of our previous paper.E2l The first ques- 
tion which we address is the observability of the four- 
fold oscillations in the density of states. These oscilla- 
tions Jaave not been seen in the experiments by Moler 
et aln; nor have they been found in recent measure- 
ments on very high quality single crystals of YBCO.cl It 
seems likely that the estimate of E a b ~ 10K-T -1 / 2 from 
the purely two-dimensional model is too high, and— the 
three-dimensionality reduces the effect significantly.ta It 
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is also possible that the orthorhombicity, which shifts 
the minima in the density of states away from the 7r/4 
directions, combined with twinning of the crystals used 
in both experiments reduces the observable anisotropy 
significantly.K3 However, even in this case, the in-plane 
anisotropy for the fields gf up to 14T used in the ex- 
periments by Wang et alB should be within the experi- 
mental resolution. A very important observation is that 
since the anisotropy in the density of states is washed out 
rapidly as the energy is increased, the in-plane anisotropy 
of the specific, heat is greatly reduced with increased 
temperature]!!^ as seen in Fig. 7 the reduction is more 
rapid for the Doppler shift density with the larger weight 
at low energies. We only consider here the possible situ- 
ation when the configuration of the vortex lattice is iden- 
tical for the field along the nodal and the anti-nodal di- 
rections; then the limiting behavior for the specific heat 
with the field along an anti-node (a = 0) and along a 
node (a = 7r/4) is easily obtained from Eqs.(|55|)-(p6|), 
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FIG. 7. Anisotropy in the specific heat between the nodal 
(dashed line) and anti-nodal (solid line) directions for different 
models. 



Therefore while the amplitude of the y/H term con- 
firms the estimates for the nodal velocities vt and va, 
and therefore for the energy scale Eh, such a term has 
not been observed for the field in the plane, and therefore 
there is no direct measurement of E a b available. How- 
ever, an estimate for this scale can be obtained from the 
scaling plots for the specific heat. 



B. Scaling functions 

It has been pointed out by Simon and LeeEl that on 
general grounds the thermodynamic coefficients of the 
nodal fermions in a magnetic field should scale with the 
variable T/Eh', consequently the experimental results 
can be interpreted as giving the form of these scaling 
functions. The scaling of the specific heat itself fol- 
lows easily from Eq.(ffl|) for the density of states, and 
the weak and the strong field limits of the scaling func- 
tion are obtained from the equations for the specific heat 
above. For the field H||c, we define Z — T/Eh and 
F C (Z) = nv f v A C(T,H)/{2TE H ); then 



F C {Z) 



dxx 2 F, 



N 



cosh 



(70) 



with Fff given by Eq. 
function follow easily: 



3). The limits for the scaling 



F C (Z) 



r 7 r 2 M 1 /3 + 77r 4 P(0)Z 2 /15, if Z < 1; 
\9C(3)Z + Z- 1 ln2, ifZ>l. 



(71) 



The numerically determined scaling function is shown in 
Fig. 8. It is remarkably similar to the scaling plot ob- 
tained from the measured specific heat in Ref. pi In that 
experiment the crossover scale, marking the transition 
from the field dominated regime, where Fc(Z) ~ const, 
to the temperature dominated regime, has been deter- 
mined to be T/y/~H ss 6.5 K/T 1 / 2 ; a very close value 
has heen obtained in a more recent experiment of Wang 
et alB. As is clearly seen from Fig. 8 the value of the 
scaling variable at the crossover depends on the struc- 
ture of the vortex state; this is easy to understand from 
Eq. ( fn] ) . The zero temperature value of the scaling func- 
tion is determined by the first moment of the Doppler 
shift distribution Mi, while the increase of Fc with the 
temperature is proportional to the weight of the distribu- 
tion at the vanishing Doppler shift, P(0). Consequently 
the crossov er value, Z c , can be expected to be propor- 
tional to \J 'Mi/ 'P(0). As the number of zeroes of the 
superfluid velocity grows, the weight in P(x) is shifted 
towards lower energies, so that Mi decreases while P(0) 
increases; these opposing trends lead to significant vari- 
ations in Z c . From Fig. 8, for the liquid and single vortex 
models the crossover occurs around Z c ~ 0.2 — 0.3; tak- 
ing this value as the experimentally determined crossover 
point, we arrive at Eh/Vh ~ 30 K-T -1 / 2 , in agreement 
with our previous estimate. Notice that the crossover 
occurs at Z c <C 1; this is simply the result of a large co- 
efficient of the Z 2 term in the low temperature expansion 
in Eq. @, 7tt 4 /15 ~ 45, while tt 2 /3 ~ 3. 
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FIG. 8. Scaling finction for the specific heat. For con- 
creteness the crossover values Z c have been defined as the 
point of a 20% increase above the high field flat region: 
F C (Z C ) = 1.2F C (0) 

Similar analysis can be carried out for the field ap- 
plied in the plane by introducing Z ab = T/E ab and 
Fab{Z ab -.a) = irv f v A C{T,H;a)/{TE ab ); the limiting 
form of the scaling functions for Z ab -C 1 can be read 
off Eqs.@-@; at Z ab > 1 we have 



F ab = 18C(3)Z ah + Z- 1 ln2. 



(72) 



The specific heat data of Refs. g|j are analyzed by mod- 
eling and subtracting the "background" contributions to 
the specific heat (phonons, Schottky anomalies etc.). To 
avoid the extensive analysis, Revaz et alu have looked 
at the difference between the specific heat with the field 
along the c-axis, and the field along the anti-nodal di- 
rection, the c/a - b difference 6C(T,H) = C{T,H) - 
C(T, H;Q). The comparison between the results of Rcf. 
f| and Refs. ||,|| has been a subject of some controversy, 
most clearly stated in Ref. [|. It has been argued already 
by the present authors and Carbotte that the experimejjp 
tal results from these groups are in fact in agreement^, 
and here we elaborate further on the sources of the appar- 
ent differences. We interpret SC(T, H) as a pure vortex 
quantity, ignoring the possible elastic contribution of the 
vortex lattice and the possible field dependence of the 
anisotropy i] — E ab /En- The issues raised in Ref. || in- 
clude the temperature dependence of SC(T, H)/T in the 
regime where C (T, H ) /T is essentially insensitive to tem- 
perature, and a form of the scaling function for SC which 
is quite different from that of C(T,H). 

As is clear from Fig. 8 for the field H||c the ratio 
C(T,H)/T does not depend strongly on the tempera- 
ture for T < Th ~ 0.1 — 0.25-Eff, reflecting the energy 
independence of the density of states for ui <C Eh- For 
the field applied in the plane along the anti-nodal direc- 
tion the physics is very similar, up to rescaling of the 
energies, which means that the density of states is only 



constant for u < E ab < E H , and therefore C(T, H; 0)/T 
is T-dependent above T ab ~ (0.1 - 0.25)E ab . The dif- 
ference, SC/T, becomes temperature dependent at the 
lower of the two crossovers, which is at T ~ 0.1E ab , and 
for T ab < T < Th it varies with temperature even though 
C(T, H)/T is approximately constant. 

It is easy to understand the difference in the scaling 
behavior between SC and C(T,H). Taking the ratio 
Eh I E a b = 4 we plot the corresponding scaling functions 
in Fig. 9. Even in the regime where the scaling functions 
Fc is nearly constant F$c is decreasing continuously. We 
therefore believe that there is no contradiction between 
the results of Ref. 0, and Refs. gjjj. Both the temperature 
dependence of SC/T and the difference in the behavior 
of the scaling function reflect the smaller Doppler energy 
scale in the plane, E ab <C Eh, and the results of these 
experiments are, in fact, quite consistent. Notice that 
the crossover in SC is much wider than that in C(T,H) 
because of two energy scales contributing to it: it ex- 
tends over a decade in the scaling variable. Note that in 
Fig. 9 we have evaluated the specific heat with the field 
H||c in the single vortex approximation, while the spe- 
cific heat for the field along the anti-node in the plane has 
been evaluated for the liquid I distribution, to model the 
expected difference in the degree of order in the vortex 
lattice. Both quantities have been evaluated with the 
single-vortex distribution in a prior publication,E3 and 
there are no qualitative differences between the two cases. 




FIG. 9. Scaling functions for the specific heat with the field 
applied along the c-axis, C(T,H) and the c/a — b difference 
SC = C(T, H) - C(T, H- 0). The former has been evaluated 
for the single vortex distribution of the Doppler shift, the 
latter for the liquid I model as explained in the text. The 
behavior remains essentially unmodified for other forms of 
the distribution. Inset: difference between the nodal (top) 
and anti-nodal directions disappears on the scale of the larger 
graph. 

To further quantify these considerations we note that 
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even though the crossover to the temperature insensi- 
tive SC has not been found in Ref. ||, the data suggest 
that it is close to T/\fH ~ 0.5 K T" 1 / 2 , which is the 
lowest value of the scaling variable reached in the paper 
(experimental measurements are limited to the temper- 
atures above ~ 1.5K, since at lower temperatures the 
non- vortex contributions to C(T,H) become dominant). 
Taking this number as a crossover value of T/E ao , we 
estimate E ab /\fH - 3 - 4.5 K T" 1 / 2 ; 2-3 times smaller 
than the estimate from the London model. 

If the value of E ab is low, it is not surprising that the 
in-plane anisotropy between the nodal and the anti-nodal 
directions has not been found in the experiments of Ref. 
|: at 14 T, E ab w 11 - 17 K, and even at the lowest tem- 
perature where the measurements of Ref. [s] have been 
made T/Eab > 0.09 — 0.15. Then the anisotropy in the 
density of states is significantly reduced from the T = 
value, see the inset of Fig. 9. On the other hand, the 
data of Ref. || for the field in the plane yield (after the 
subtraction of the Schottky anomaly) a crossover temper- 
ature between the field dominated and the temperature 
dominated regimes close to T cr ss 2 K per T 1 / 2 . If this 
value is taken as corresponding to the the crossover in 
T/E ab , it implies a large value of E ao /\fH > 10 — 20 K 
T -1 / 2 . In that case the absence of the anisotropy can 
only be explained (somewhat unsatisfactorily) by an ap- 
peal to the three-dimensionalityli3-at-.a combination of 
the orthorhombicity and twinningJlj'E3 Since part of the 
experimental difficulty stems from the smallness of the 
>/H term with the field in the plane, the analysis for 
that geometry typically involves^assurning a field depen- 
dent contribution of that format! However, if E ab is 
small, the field dependence of the specific heat is modi- 
fied by the Zeeman splitting, and this splitting has to be 
taken into account in the analysis. 

C. Zeeman splitting 

If the energy scale for the in-plane Doppler shift is 
indeed much smaller than the naive estimate from the 
London model, E ab /y/~H ~ 3-4.5 K T" 1 / 2 , the Zeeman 
splitting has a significant effect on the specific heat with 
the field applied along a node in the experimentally rele- 
vant range. The specific heat no longer obeys the scaling 
properties discussed above; for the field along a node the 
contribution of the Doppler- "inactive" nodes is given by 

C z (T,H) = -^-F z (f^-) (73) 

F z {x)=x t 2 cosh" ~ 2 -dt+ / i 3 cosh" 2 -eft, (74) 
Jo 2 J x 2 

in agreement with Ref. |6(], and therefore scales with H 
rather than y/~H. 

As the in-plane anisotropy in the density of states has a 
maximum for u> = [J,H, the anisotropy in the specific heat 



also goes through a maximum; we expect approximately 
T m ax(H) oc H. We consider here two different cases for 
a = E ao I \J~H: a large value corresponding to our original 
estimate a — 10K T" 1 / 2 , and a small value implied by the 
experiment, a = 4K T -1 / 2 , and evaluate the specific heat 
for the liquid I distribution. The main panel of Fig. 10 
shows the scaling plot for the in-plane anisotropy in the 
specific heat, C ams (T,H) = C(T,H;0) - C(T,H;n/4) 
at H = 10T, so that the values for the two cases are 
E a i, ~ 32K, and E ab ~ UK. While the anisotropy is 
severely reduced at T = Q, it becomes close to the val- 
ues estimated without accounting for the Zeeman shift 
at the temperatures T ~ 2.2K and T ~ 1.8K for the 
two cases respectively, and therefore the anisotropy in 
the experimentally relevant regime is not modified sig- 
nificantly. Nevertheless, if the absolute magnitude of the 
anisotropic term is small, and its field dependence has to 
be modeled in the analysis of the experimental resultsQ, 
it is important to note that, as is clear from the inset 
of Fig. 10, the field dependence of the anisotropy is not 
simply proportional to \/H, but flattens and decreases 
at high fields. The deviations are especially important 
for small a, since then the maximum of the anisotropy is 
reached at H ~ 10T for T = 1.5K, well within the exper- 
imental range. We analyze this scenario in more detail 
in Fig. 11, which demonstrates that if the coefficient a is 
small, the maximum in the anisotropy can be observed 
at low temperatures, but moves out of the easily acces- 
sible range, to H > 20 T, at higher T. In comparison, if 
a r*> 10, the maximum lies at high fields for all relevant 
T. 




FIG. 10. Main panel: the anisotropy in the in-plane spe- 
cific heat plotted in the scaling form at H = 10 Tesla for 
the coefficient a = 4, 10 K T -1 / 2 , and without accounting for 
the Zeeman shift (dashed line). Inset: The specific heat at 
T = 1.5K, which is close to the lowest experimentally accessi- 
ble temperature, for the same two values of a with (solid line) 
and without (dashed line) accounting for the Zeeman shift. 
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FIG. 11. Anisotropy o fthe in-plane specific heat for a = 4 
K T -1 / 2 as a function of the field at different temperatures. 

Consequently, in the search for the experimental verifi- 
cation of the anisotropy in the specific heat, it cannot be 
assumed that the anisotropy increases as \[H\ if the en- 
ergy scale for the in-plane Doppler shift is small, the Zee- 
man splitting modifies the field dependence of the specific 
heat. For the small value a = 4K T -1 / 2 the maximum 
anisotropy, reached in the fields of the order of 10-15 
Tesla at T = 1.5-3 K, is of the order of 0-5rfjhft7s; based 
on the available experimental values for 7 s 0tl'El between 
0.064 and 0.21 mJ mol -1 K -3 , the maximal anisotropy 
ranges between 0.032 and 0.19 mJ mol -1 K _1 ; it is sig- 
nificantly larger for larger values of E a b / V~H. 

Recently, Wang et al. attempted to observe the an- 
gular oscillations we have predicted in the in-plane spe- 
cific heatB They did not, however, find appreciable dif- 
ference between two measurements with field applied in 
the nodal and anti-nodal directions. Several reasons may 
have contributed: first of all, the YBCO sample used in 
their experiment is twinned. Twinning, combined with 
the orthorlpmbicity of YBCO is expected to reduce the 
anisotropy.E3E3 Here we point out another possible rea- 
son for the difficulty in extracting the difference between 
the two directions from the data: the field dependence 
of the anisotropic term is not simply given by yH, as it 
would be in the absence of the Zeeman term, and as as- 
sumed in Ref. 0. Instead, the anisotropy increases with 
the field up to fields of about 10-15 T, and decreases 
thereafter. Consequently, we believe that to confirm the 
predicted oscillations experimentally, it is highly desir- 
able to use an untwinned crystal, and carry out the mea- 
surements at intermediate fields (10-15 T) at the lowest 
possible temperatures, since the anisotropy in the specific 
heat is expected to be the largest in this range. 



VII. SPIN-LATTICE RELAXATION RATE 

We now turn our attention to the calculation of the 
response functions. In these calculations the local, in real 
space, physical quantities depend on the Doppler shift at 
both pairs of nodes, and consequently the averaging has 
to be carried out with the two-node probability density C 
rather than the single node distribution V . The simplest 
example of such a quantity is the average spin-lattice 
relaxation rate which we now consider. 

Since the NMR measurements on cuprates are typi- 
cally done in a magnetic field of ~ 10T the effect of the 
field on the measured signal has to be considered in the 
analysis of the data. There are at least two effects of the 
vortex state on the spin-lattice relaxation time. First, the 
Doppler shift modifies the local density of states, intro- 
ducing the local relaxation rate, which varies from point 
to point. Second, the magnetization due to the vortex 
lattice introduces inhomogeneities in the field, leading to 
the broadening of the resonance line. As a result, there 
are two possible approaches to the analysis. In a perfect 
vortex lattice there exists a one-to-one correspondence 
between the local field at a particular point in the unit 
cell of the vortex lattice, and the value of the superfluid 
velocity at that point. Assuming such a perfect lattice it 
is therefore possible to associate the local relaxation rate 
with the relaxation rate at a particular frequency in the 
resonant line. Such an approach has beep .developed the- 
oretically in the semiclassical frameworkEilE 2 ], and the re- 
sults are in qualitative agreement with the experimental 
observation that the relaxation rate and the local den- 
sity of states are larger in the regions of higher field, i.e. 
higher supervelocityu. 

On the other hand, in a disordered vortex state there 
is no unique identification of the local value of the su- 
pervelocity corresponding to a local magnetic field. It 
may therefore be useful to analyze the relaxation rate ob- 
tained by the "global" fit to the resonance line; especially 
when the linewidth remains quite narrow in frequency. It 
has been shown that the time-decay of the magnetization 
is non-exponential as it involves a convolution of many 
local relaxation rates, but that it is possible to describe 
it with an effective scatteriugxate which depends on the 
field and the temper ature.EjH Usually the analysis of 
the experimental data is done assuming a single relax- 
ation rate, and it is therefore important to understand 
its behavior in a d-wave superconductor. 

The analysis of the relaxation time, T±, can be under- 
taken either by looking at its magnitude directly, or by 
analyzing the ratio Ti/t c , where r c is the relaxation time 
at T — T c . The former approach involves modeling or 
estimating from the available data the matrix element 
for the interaction; it has been used, for example, in Ref. 
|5l| . The latter method is based on making assumptions 
about the normal state relaxation in the cuprates. We 
employ it assuming a normal metallic relaxation at T c 
with the caveat that this may be only qualitatively cor- 
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rect for underdoped compounds. 

With this assumption for the spin-i system the mag- 
netization decays as m(t) = M(t)/M(0) — exp(— t/Ti), 
where the relaxation rate in the infinitesimal field is given 
by 



t c T c 
T X T 



4-oc 



duj 



duo J 



1 

2./o 

.2 



+ OC 



dx 



N 2 {xT) 



cosh x/2, 



where Nq — m/2-Kft, is the 2D density of states in the 
normal state. 

In non-zero magnetic field, the decay of the average 
magnetization is given by 



m{t) = 4 



poo poo 

/ de 1 / de 2 £(ei,e2)exp[-t/Ti(ei,e 2 ) 
Jo Jo 



(76) 



where the position and Doppler shift dependent relax- 
ation rate is determined from 



t c T c 



1 



2iV 2 



dxN 2 {xT,e 1 ,e 2 )coshT 2 x/2. (77) 



The density of states is given by the sum of the contribu- 
tions of all nodes, \lu + {i^v fv&) with (ej = ±e l7 ±e2), 
which yields 



N(u,e 1 ,e 2 ) = r— ^ 

( 2w, 



(78) 



if lu > max(ei, e 2 )] 
x < lu + max(ei, e 2 ), if min(ei, e 2 ) < lu < max(ei, €2); 
[f-2+ei, if lu < min(ei, e 2 ). 

Here again, without loss of generality, we set u>, ei, £2 > 0. 

The decay.|_of the magnetization is therefore non- 
exponentia£jl£j, as is shown in Fig. 12 for a liquid dis- 
tribution. In a magnetic field the density of states is 
enhanced globally, and therefore the stronger the field, 
the faster the decay of m(t). In the regime Eh 3> T the 
density of states is significantly enhanced over a large 
part of the unit cell of the vortex lattice. This high den- 
sity of states yields a fast relaxation rate responsible for 
the initial decrease in the magnetization. The long-time 
decay of m(t) is determined by the slowest relaxation 
rates, which occur in the regions where the superfiuid 
velocity is small and the density of states is largely de- 
termined by the temperature. The two regimes are seen 
in Fig. 12 : the large-i tail of lnm(i) is affected by the 
field much more weakly than the short-time decay. For 
all the values of the field it is possible to fit the time 
dependence by an exponential, although clearly the re- 
laxation rate obtained from such a fit differs significantly 
from the zero-field rate. We have addressed the fit of 
the magnetisation at different time scales in a previous 
publicationJla'Ea 
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FIG. 12. Magnetization decay at a fixed temperature 
for different values of the magnetic field. We have used 
2No~KVfVA = kfVA ~ 2Ao (pure d-wave), and have set 
Ao = 2.14T C . m(t) has been evaluated for the liquid I model. 

An important comment concerns the scaling of the 
magnetization. First of all, due to scaling properties 
of the density of states, the magnetization decay due to 
spin-lattice relaxation satisfies 



m(t)=F m (tHTf( 



T 
H 



(79) 



where the functions F m and / can be obtained from the 
general expression Eq.(f76|). Moreover, when Eh/T ^> 1 
the density of states and therefore the function / are 
nearly constant. Two conclusions follow immediately. 
First, at a fixed ratio T/\/~H the magnetization depends 
only on the single variable tTH. Second, at low tem- 
peratures m(tTH) is independent of the ratio T/yH 
at short time scales, when the relaxation rate is domi- 
nated by the field-induced density of states rather than 
the temperature-driven density of states. The collapse of 
the low-T data on a Wfile curve as a function of tT has 
been found previouslyoO, however we are not aware of an 
experimental check of such scaling at different fields. In 
Fig. 13 this behavior is clearly seen. Deviations from the 
scaling form are noticeable already for Eh/T ~ 4, how- 
ever even in this regime the curves for different Eh and T, 
but with the same ratio Eh/T, coincide. The scaling is 
always obeyed at short time scales, where the time decay 
of the magnetization is determined by the fast relaxation 
rates in the regions with a large Doppler shift. On the 
other hand, at long time scales the time-dependence of 
m(t) is determined by the slowest relaxation rates, in the 
regions where the Doppler shift vanishes, and therefore 
there are always deviations from the scaling with tTH. 
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FIG. 13. Magnetization as a function 
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The ratio T/Eh is identical for the 



Since the long time scale decay is determined by the 
measure of the points with small Doppler shift, it de- 
pends crucially on the probability density £(ei,£2). In 
particular, there is a dramatic difference between the sin- 
gle vortex picture, where e 2 + e 2 > Ejj, and the lattice 
or liquid states, where this restriction is lifted: magne- 
tization decays much faster in the single vortex picture, 
as can be seen in Fig. 14. Notice that for the very early 
times, when the magnetization decay is determined by 
the regions with the highest Doppler shift, the two dis- 
tributions give the same result. Therefore the effective 
relaxation rate obtained from the exponential fit depends 
not only on the field but also on the structure of the vor- 
tex state. 
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FIG. 14. The difference between the liquid (liquid I) and 
the single vortex models in a strong (main panel) and weak 
(inset) magnetic field. Parameters are the same as in Fig. 12. 

The difference in the behavior for the two types of the 



vortex state can be understood from the analysis of the 
magnetization. Noticing that £(£1,62) = £i(ef + due 
to symmetry ( the probability density should be even 
in both £1 and £2, and should be symmetric under the 
interchange £1 «-> £2), and introducing polar coordinates 
x = (e 2 +€%)/E%, tan</> = ei/e2, we arrive at 

/>OC 

m(t) m / dxC^e'^IiPx/A), (80) 

Jxq 

r 

I(z) = / exp(— zsm(p)d4>, (81) 
Jo 

where x w T 2 /E 2 H , and /3 = tTE 2 H /t c T c AI. This form 
shows explicitly that there is an approximate scaling with 
the variable 0, and that the scaling is obeyed better the 
smaller the ratio T/Eh- 

Due to the exponential the integral over x is cut off at 
j3x 1, so that we only need to evaluate I(z) for z ~ 1. 
In that case all angles <fi contribute to the integral, and 
I(z) w 7rexp(— bz) with b ~ 1 (in contrast, I(z) ~ 2/z 
for z ^> 1), leading to 



m(t) w 7T / dxCi(x)e 



-Pax/4 



(82) 



with s = 1 + a ~ 1. Consequently, the long time scale 
limits (/3 > 1) for the single vortex and liquid I distribu- 
tion respectively 
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(83) 



(84) 



As a result, the long time decay for the liquid regime is 
governed by the relaxation rate close to the (3xo oc T 3 
behavior expected for H = 0, while for the single vortex 
model the relaxation rate is proportional to f3 oc TH . 

In reality however the decay of m(t) at long time scales 
is usually not measured, and at intermediate times the 
detailed analysis of the time dependence of the magne- 
tization taking into account the non-exponential form of 
m(t) is complex. It is possible to define an effective re- 
laxation rate, however, the weight of the components of 
the magnetization with fast and slow decay is different 
for different definitions, and the resulting effective relax- 
ation rate is different, as we now illustrate. One possible 
approach is to define the effective rate as 
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eff 



= 4 



dei / de 2 £(ei,£2) 
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(85) 



Unlike the average for the magnetization, Eq.(|7q), which 
has the largest contribution from the slowest relaxation 
rates, in this averaging procedure the weight of short 
relaxation rates is high, and a cutoff of the energy integral 
near the core is required. To leading order in In E q /Eh 
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the relaxation rate in the field dominated regime is then 
given by 
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7T + 2 1 



eff 
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(86) 



The relaxation rate given by this expression is expected 
to overestimate the rate of the decay of the magnetiza- 
tion. Fast relaxation occurs near the cores, where the 
effective field is higher, and therefore in the component 
of the signaL-away from the original position of the res- 
onance line.EJ Alternatively, we can define the average 
relaxation time 



*ff 



m(t)dt 



(87) 
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This average has a large contribution from slow relax- 
ation rates, and we expect the effective rate to be under- 
estimates since over experimentally relevant time scales 
the slowest rates do not contribute to the magnetization 
decay appreciably. Indeed, for the cases of the single 
vortex and the liquid distributions we obtain in the field 
dominated regime 
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(88) 
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The coefficient in the last expression is significantly 
smaller than the expression given by Eq.([S6|). We can 
now compare this expression with the result of Ref. |^, 
where it was found that T c T c jT\T w 0.2 at H = 11 T 
at low T. From our estimate of Eh it follows that at 
this field E H ~ 100K. Taking E ~ 1500 K, we find this 
ratio for the average rate to be be 0.35, while the average 
relaxation time procedure yields the values of 0.06 and 
0.005 (at T ~ 5 K) respectively. The experimental value 
is between the two estimates, as expected. 



VIII. IMPURITY SCATTERING 

A. The self-energy 

In the presence of impurity scattering the frequency 
is renormalized according to ui = w — £(£). The self 
energy £(£>) depends on the momentum integral of the 
Green's function, and therefore on the Doppler shifts at 
all nodes. Consequently, in all the calculations involving 
the impurities, the local quantities depend on both ei 
and e 2 . 

Here we consider the impurity scattering in the unitar- 
ity limit. The strategy for the calculation is as follows. 



The self energy is given by Eq.(g); to evaluate this expres- 
sion in a field we have to introduce the Doppler shift in 
the Green's function as before, and solve for the self en- 
ergy self-consistently at each node. In other words, there 
is a distinct Doppler shift at each node, and the self- 
consistency requires that the scattering to other nodes 
with their respective Doppler shifts be taken into account 
self-consistently. Therefore we can write 



£(w, ei,e 2 ) = -ni[Go(w,ei,e 2 )]" 



where, from Eq. (@). 
G (w,ei,e 2 ) = - 
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(92) 
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We now focus our attention on the cases of weak 
u> <C Eh <C 7, and strong ui <C 7 <§C Eh fields at low 
temperatures. Setting u = it is clear immediately that 
the real part of the momentum integral of the Green's 
function at each node in Eq. ( |9l| ) is odd in the Doppler 
shift, e„, and therefore vanishes upon summation. t2l As a 
result, the self energy has only the imaginary part given 
by Eq.® with 
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In the strong field limit, cj 2 <C £i,e 2 , we obtain for the 
the density of states 
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as expected, see Eq.(^8 
this regime is given by 
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(96) 



In the weak field impurity dominated regime, a> 2 3> ei, e 2 , 
on the other hand, mv the field-induced change in the 
density of states is quadratic in the Doppler shift and is 
given by 



<HV(ei,e 2 ) 
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(97) 
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where 7 is the zero en ergy scattering rate which has 



been defined in Section [I|. Then to the leading order 



in In Eq I En the average change in the density of states 
is given by 



ZTT^-fVfVA Eh H 



(98) 



for the single vortex and the liquid distributions. For 
the Gaussian model the change in the density of 
states is smaller by a logarithm of a large number, 
SN s (0,H)/SN G (0,H) = 2\nE /E H . This behavior is 
illustrated in Fig. 15. Notice that the three distributions 
yield different high-field slopes, corresponding to the dif- 
ferent values of the moment Mi . 




[Telia 1 ' 2 ] 

FIG. 15. Residual density of states as a function of the 
field. The impurity scattering is taken in the unitarity limit, 
with (runVfVA) 1 ^ 2 = 20 K. The parameters are Eq = 1500K, 
E H /V~H = 30 K T- 1/2 . 

Theas-. results are in agreement with the previous 
worklijUj. The low energy scattering rate, 7, provides 
the new energy scale in addition to the average Doppler 
shift Eh and the temperature T. At low temperatures, 
the competition between Eh and 7 determines the be- 
havior of the density of states, and in the field dominated 
regime, Eh S> 7, the density of states strongly depends 
on the probability density of the Doppler shift, as it does 
in the pure limit. The dependence of the self energy on 
the magnetic field is crucially important for the analysis 
of the transport properties in the vortex state, and we 
will discuss these issues in detail in a separate paper. 



IX. CONCLUSIONS 

In this paper we have discussed the semiclassical ap- 
proach to the vortex state of unconventional supercon- 
ductors and have applied it to the analysis of the ther- 
modynamic properties of a two dimensional <i-wave su- 
perconductor, which we take as a model for the high-T c 



cuprates at low energies. Our main point is that within 
the semiclassical approach the dependence of the mea- 
sured quantities on the magnetic field is sensitive to the 
structure of the vortex state and the distribution of the 
supercurrents. This is shown in an approach which in- 
volves introducing the Doppler shift due to circulating 
supercurrents into the quasiparticle dispersion, and com- 
puting the physically measured magnetic field dependent 
quantities as a spatial average of their local values in the 
vortex state. The major step which has enabled us to 
move beyond the standard single vortex description is 
the rewriting of the spatial average in terms of the aver- 
age over a probability density of the Doppler shift at a 
particular node or at a pair of nodes. We have analyti- 
cally computed these probability densities for the single 
vortex picture, for model liquid distributions, and for a 
non-physical, albeit often used, completely random dis- 
tribution. We have argued that this approach is easily 
applicable to any given distribution of vortices, and that 
the single vortex and the liquid models typically give the 
upper and the lower limits of the field dependence since 
they over- and under- estimate respectively the number 
of the points in the vortex lattice unit cell where the 
Doppler shift dominates the physical picture. 

We have applied this approach to the analysis of the 
electronic specific heat in the vortex state and to the de- 
scription of the spin-lattice relaxation of the NMR mag- 
netization. In the former case the specific heat depends 
on the single-node probability distribution. The values 
for the Fermi velocity at the nodes, as well as of the slope 
of the superconducting gap, determined from such an 
analysis are consistent with the values inferred from other 
experimental measurements, and the values directly de- 
termined from the photocmission. Moreover, noticing 
that the magnitude of the y/~H term is larger for the more 
ordered vortex state, has allowed us to reduce the discrep- 
ancy between the results for the gap slope obtained by 
different experimental groups. We have also emphasized 
that the difference in the form of the scaling function 
obtained by these groups is naturally explained as a con- 
sequence of the smaller Doppler shift energy scale for the 
field applied in the plane; this work confirms our eaxlier 
assessment on the basis of the single vortex pictureS 

Since the analysis of the scaling plots allowed us to es- 
timate the energy scale of the Doppler shift for the field in 
the plane, and since this energy scale is smalles than the 
London model estimate of our previous worklla, we have 
investigated here whether the anisotropy in the specific 
heat between the experimental arrangements with the 
field applied along a node and between the two nodes is 
observable. We have paid special attention to the effect of 
the Zeeman splitting, which becomes more important for 
smaller in-plane scale E a f,. We have found that, while the 
zero temperature anisotropy is significantly reduced com- 
pared to the case of no Zeeman splitting, as predictecO, 
the anisotropy does not decrease with the temperature 
and in the experimentally relevant temperature range the 
magnitude of the anisotropy is weakly affected by the in- 
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elusion of the Zeeman splitting. The field dependence 
of the anisotropic specific heat may however be modified 
quite significantly, and this change has to be taken into 
account when analysing the experimental data. 

We have considered the spin-lattice relaxation rate as 
an example of a response function which depends on the 
probability distribution at two nodes; in contrast to the 
specific heat the contributions of the nodes are not sim- 
ply additive. It has been known that the magnetization 
decay is non-exppaential due to a distribution of the local 
relaxation times .Ej We have shown here that the effective 
relaxation rate obtained from a fit to an exponential at 
short or long time scales is different, and at long times 
depends crucially on the structure of the vortex state. 
We have predicted a scaling form of the magnetization, 
and discussed the existing evidence for such a scaling; 
of course more experimental work checking this predic- 
tion would be highly desirable. We have also introduced 
an effective relaxation rate obtained from a global fit to 
the magnetization, and found that it agrees qualitatively 
with the available experimental results. 

In general, the structure of the vortex state may change 
quite dramatically as the temperature and the applied 
field are varied; one example of such a change is the melt- 
ing of the vortex lattice. In such a situation we expect a 
change in the spatially averaged thermodynamic quanti- 
ties measured in experiment which reflects the transition 
from one type of distribution to another. Moreover, as 
the degree of ordering in the vortex lattice depends on 
the history of the sample, the measured field dependence 
varies accordingly. For example, the coefficient of the 
\fH term in the electronic specific heat should, in gen- 
eral, depend on whether the sample has been cooled in 
an applied field or in zero field: in the latter case the vor- 
tex state is more disordered. Whether these effects are 
observable experimentally depends crucially on the qual- 
ity of the sample since the changes may be rather small, 
nevertheless, in clean untwinned samples they may be 
measurable. 

Finally, to illustrate that in the presence of impuri- 
ties the two-node probability density is always required 
we have analysed the density of states in the impurity- 
dominated regime for different structures of the vortex 
state. This part of our work will be developed further in 
the analysis of the transport properties, which warrants 
a separate paper. 
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